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\Q \ Abstract. Based on detailed 2D radiation hydrodynamics (RHD) sim- 

ulations, we have investigated the effects of photospheric temperature 
inhomogeneities induced by convection on spectroscopic determinations 
of the lithium abundance. Computations have been performed both for 
t~^. \ the solar case and for a metal-poor dwarf. NLTE effects are taken into 

<^> • account, using a five-level atomic model for Li I. Comparisons are pre- 

sented with traditional ID models having the same effective temperature 
and gravity. The net result is that, while LTE results differ dramatically 
between ID and 2D models, especially in the metal-poor case, this does 
not remain true when NLTE effects are included: 1D/2D differences in 
the inferred NLTE Li abundance are always well below 0.1 dex. The 
present computations still assume LTE in the continuum. New computa- 
tions removing this assumption are planned for the near future. 
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1. Introduction 

It would be an offense to the audience here to pretend to explain why it is 
important to determine accurately the abundance of 6 Li and 7 Li in the oldest 
stars. In this respect, we have nothing to add to the exposition by Erangois 
Spite (this volume). We will right away mention the two major problems which 
may cast doubts on our real knowledge of the actual initial abundance of Li in 
the oldest stars, and consequently in the primordial matter, (i) While standard 
models of the internal structure of metal-poor dwarfs do not deplete 7 Li, more 
sophisticated models including rotationally induced mixing (Pinsonneault et al. 
1992) have predicted that the measured abundance in the photosphere is 5 to 10 
times less than the initial abundance representative of Big Bang material, (ii) On 
top of that, Kurucz (1995) claimed that the hot and cold convective structures 
produce large effects in metal-poor stellar photospheres, where the convection 
zone reaches the line formation layers. The claimed effect is an overionization of 
Li by a factor of 10, leading to an underestimation of the abundance of Li when 
derived from the resonance line of Li I (A 670.8 nm) in the usual way. 
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If these two statements are correct, the true abundance of Li in primordial 
matter is 50 to 100 times higher than the value derived from ID, LTE models of 
halo subdwarfs so far. The first factor of 5 to 10 has been discussed in a previous 
paper by Ryan (this symposium), and shown to be likely much smaller, of the 
order of 1 to 1.4. We shall not come back to this point, which we consider as 
very well treated. 

Before this symposium, a single paper (Asplund et al. 1999) has dealt with 
the question of the other factor of 10 claimed by Kurucz (1995), whose arguments 
were based on a simplified two-column model. In contrast, the work by Asplund 
et al. relies on realistic 3D hydro dynamical models, similar to the simulations of 
the solar granulation (Stein & Nordlund 1998), but with parameters appropriate 
for two metal-poor stars: HD 140283 and HD 84937, both subgiants. The 
computation of the lithium resonance line was made under the assumption of 
LTE, and the correction to be applied to the Li abundance derived from standard 
ID models was found to be large, of the order of -0.2 to -0.35 dex. Note that these 
corrections have the opposite sign as Kurucz's prediction! However, Kiselman 
(1997, 1998) had shown, in the solar case, that NLTE and LTE computations 
lead to significantly different values of equivalent widths of the Li I A 670.8 nm 
line over hot and cold structures (see Fig. 3 of his 1997 paper, top panel). 

For this reason, we decided to undertake NLTE radiation hydrodynamics 
computations for the case of a metal-poor star, and we report here on the re- 
sults of this investigation. In the next section we recall former work related to 
simulations of the solar granulation, a useful benchmark for checking the theory, 
but not directly applicable to metal-poor stars. In section 3 we describe the 
assumptions underlying the construction of the 2D RHD models used for the 
spatially resolved computation of the lithium resonance line. Section 4 gives the 
description of the NLTE treatment of the Li atom, and section 5 summarizes our 
results and compares them to those presented by M. Asplund (this symposium). 
Finally, our conclusions are listed in section 6. 

2. Former work at solar metallicity 

While there is only one paper dealing with multidimensional atmospheres for 
metal-poor stars (cited above), there are several studies for the solar case, aimed 
at understanding the variation of the continuous radiation intensity (granula- 
tion), and the behavior of spectral lines across the solar granulation pattern. 
Several of these works use snapshots from 3D simulations by Stein & Nord- 
lund (1998), such as Kiselman (1997,1998) and Uitenbroek (1998). Gadun & 
Pavlenko (1997) use their own 2D simulations. 

Of particular interest for us are the papers dealing with the combined effects 
of multidimensional structures and NLTE (Kiselman and Uitenbroek) . It is clear 
from Kiselman (1997) that the NLTE behavior of the Li I A 670.8 nm resonance 
line is drastically different from its LTE behavior, in 2D as well as in 3D models. 
The result which is the most relevant for us is the difference of 30 per cent on the 
predicted mean equivalent width (W) of the line, leading to a similar change in 
the derived lithium abundance (NLTE/LTE abundance correction +0.15 dex). 
But another interesting difference is the reverse behavior of the equivalent width 
W as a function of surface continuum brightness I c . While LTE computations 
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Metal-poor Sun ([M/H]=-2): Model OT57G44N7, t= 34140 s (TEMPERATURE) m<p) iog 1D w 




x [km] 

Figure 1. Snapshot from the metal-poor Sun simulation, showing 
the velocity field (arrows) and temperature structure (black contour 
lines). White lines are curves of iso-optical depth. The origin of the 
geometrical scale (left) corresponds to tr oss = 1; scales at right refer to 
average gas pressure P and average optical depth tr oss . This simulation 
was done on a 210x106 grid (tick marks at right and top), with a vertical 
extent of about 1 500 km (~ 7H P ), and a horizontal period of 5 250 km 
(upper and lower boundaries 'open', lateral boundaries periodic). 

result in a strongly positive slope in the W versus I c diagram (with a large 
scatter around the mean relation) , NLTE computations show a slightly negative 
slope and a much tighter (anti-)correlation between I c and W. This reflects the 
fact that the population of the Li I levels is much more controlled by the local 
temperature in LTE than in NLTE, where, for weak lines, the photoionization 
rates play the dominant role. So, even if the general conclusion of the above 
mentioned papers is that, in the solar case, the abundance determination of Li 
is not strongly affected by the combined effects of temperature fluctuations and 
NLTE, in comparison to what is obtained with classical ID models having the 
same effective temperature and gravity, it appears unsafe to compute abundances 
from multidimensional stellar atmospheres based on the assumption of LTE line 
formation. 

3. 2D radiation hydrodynamics models 

Our LTE and NLTE computations have been performed on the basis of several 
snapshots from a 2D numerical simulation of convection in a stellar envelope 
having the same effective temperature and gravity as the Sun, but a 100 times 
smaller metallicity. Basically, the time dependent equations of hydrodynamics 
are solved for a compressible fluid, with an energy equation including 3 terms: 
turbulent and shock dissipation of kinetic energy, diffusive transport of heat, and 
radiative energy exchange. The main limitation of the code is the restriction of 
the flow to two spatial dimensions. Magnetic fields and rotation are ignored. 

Apart from these simplifications, as much realistic physics as possible is 
included. The equation of state accounts for ionization of hydrogen and helium 
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as well as H2 molecule formation, opacities have been adapted from Kurucz's 
ATLAS code and include line absorption. For the computation of the radiative 
energy balance, we employ a multi-dimensional, non-local, frequency-dependent 
radiative transfer scheme, actually solving the transfer equation along 26880 
independent rays of various inclinations, using an efficient modified Feautrier 
method (Feautrier 1964). At the bottom boundary, inflowing matter has a given 
specific entropy, which is adjusted to produce the prescribed effective tempera- 
ture of the atmosphere. Energy dissipation on small scales is roughly modeled by 
introducing a subgrid scale eddy viscosity, depending on the grid resolution and 
local velocity gradients in the usual way. Details of the employed hydrodynamics 
code can be found in Ludwig et al. (1994) and Freytag et al. (1996). 

Fig. 1 shows a sample snapshot from our metal-poor Sun simulation. Note 
the complex velocity pattern and the occurrence of very strong temperature 
gradients. The relevant region for the formation of the Li I line is the r = 0.1 
contour line. 



4. NLTE computation of the Li I resonance line 

The computation of the Li I spectrum is greatly simplified by the fact that 
all lines of Li I are weak, Li being a trace element. So the radiation field in 
the line is, to first approximation, the same as the continuous radiation field, a 
single iteration being sufficient for taking care of the small perturbation of the 
monochromatic radiation field brought about by the line. 

We have, as a first step, approximated the Li I atomic configuration by a 
five level atom, exactly as done by Uitenbroek (1998). This leaves six permitted 
bound-bound transitions, and five photoionization rates needing the computa- 
tion of the continuous radiation field at frequencies above the threshold, until the 
contribution of the product of photoionization cross-section and mean intensity 
of the radiation becomes negligible. Because in the UV the contribution of lines 
to the opacity is important, we have used Kurucz's Opacity Distribution Func- 
tions (ATLAS 9) for the relevant metallicity. This multiplies the computation 
time by 12, as each opacity bin is subdivided into 12 subintervals. So, for each 
snapshot, the transfer equations must be solved for about 12 x 120 wavelengths, 
along 26880 different rays (note that these extensive computations are done only 
after the actual hydro dynamical simulation for a few selected snapshots). After 
this, it is possible to compute all the coefficients in the equations of statistical 
equilibrium (see Mihalas 1970, p. 144). Once the departure coefficients bi are 
evaluated, the line can be computed using the source function: 



Hi ~i~ H c fc>\ + K c 

where: 



2hvl 1 



&l - Zij - — 2~ 



c 2 (pi/bj) exp(huij /kT) — 1 

is the line source function. The departure coefficient for the lower level is bi and 
the one for the upper level bj. The other notations are standard. Subscript "c" 
stands for continuum, and "1" for "line" . Note that in NLTE the expression of 
the partition function is modified and becomes: 
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U = ]T hg t exp{- Xl /kT) 

i 

where the gi and Xi are the statistical weight and the excitation energy of level 
i, respectively . 

5. Results and discussion 

We have first tested our program on a Kurucz ATLAS 9 ID solar model to see 
whether the hi had the expected behavior, already computed by Carlsson et al. 
(1994). Fig. 2 shows the depth-dependence of the first 3 departure coefficients, 
applying to levels 2s, 2p and 3s, respectively. 
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Figure 2. Variation of the first three departure coefficients with 
depth (actually temperature in this graph), based on an ID ATLAS 9 
model of the solar atmosphere (mixing-length parameter a = 0.5). 
Circles refer to the ground level (b\), triangles to the 2p level (62)) an d 
crosses to the 3s level (63). As expected, all three coefficients become 
very close to 1 in the deep layers, and reach values below 0.5 in surface 
layers, as found by Carlsson et al. (1994). The adopted logarithmic 
abundance of Li is 2.2, on the scale log(nH)= 12.0. 

Next, we have computed the equivalent width of the Li I resonance line 
for the Kurucz ID solar model and for two 2D solar snapshots, still for the 
logarithmic Li abundance 2.2. Fig. 3 shows the variation of the equivalent width 
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Figure 3. Top: LTE (circles) and NLTE (crosses) equivalent widths 
of the Li I resonance line as a function of horizontal position x, for 
a representative snapshot from a 2D simulation of solar granulation. 
The thin line gives the continuum intensity at A 670.8 nm in arbitrary 
units. Bottom: Same equivalent widths as a function of continuum 
intensity ISURF. All data derived from vertical rays (^=1). 



W (for the intensity normal to the surface) of the Li I 670.8 nm line over the 
simulated granulation pattern, both as a function of the horizontal position x 
(top) and as a function of the continuum intensity I c (bottom) for one particular 
snapshot. Note the wide variation of W computed in LTE, compared to the much 
more limited excursion of W computed in NLTE. The mean equivalent widths 
for the flux spectrum integrated over the full length of the sample are given in 
Table [l] for the two solar snapshots and for the ID reference model having the 
same effective temperature, gravity and (solar) metallicity. In each case, the line 
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Table 1. LTE and NLTE mean equivalent width (W) [mA] of the 
Li I A 670.8 line for 2 different snapshots from a 2D hydrodynamical 
simulation of solar surface convection, obtained from the horizontally 
averaged flux spectrum. For comparison, the equivalent widths result- 
ing for a ID ATLAS9 reference model, computed with the same line 
formation code, are given in the first row. The results from Kieslman 
(1997) given in the last rows have been rescaled to the same Li abun- 
dance of 2.2; they refer to intensity (// = 1). 



ID model or RHD snapshot 


(W) (LTE) 


(W) (NLTE) 


ID Kurucz (a = 0.5) 


43 


32 


2D, L71D09-605 


59 


35 


2D, L71D09-625 


71 


36 


Kiselman ID (OSMARCS) 


44 


42 


Kiselman 3D 


52 


37 



is computed in LTE as well as in NLTE. We note that, in NLTE, the results for 
the 2D snapshots do not differ significantly from the ID case. 

Finally, we have carried out a similar procedure for our metal-poor stellar 
example, again computing the equivalent width of the Li I resonance line for a 
Kurucz ID model and for five snapshots from our 2D hydrodynamical simulation; 
as before, a Li abundance of 2.2 was adopted. Fig. 4 shows the variation of the 
equivalent width W of the Li I resonance line over the stellar granulation pattern 
for the typical snapshot displayed in Fig. 1. The difference between LTE and 
NLTE is even more pronounced than in the solar case, the NLTE correlation 
between W and I c being much tighter and of opposite sign compared to LTE. 

Table 2. Same as Table [l], but for 5 different snapshots from a 
2D hydrodynamical simulation of surface convection in a metal-poor 
star ([M/H] = —2) with solar effective temperature and gravity. For 
comparison, the results for the corresponding ID ATLAS9 model, com- 
puted with the same line formation code, are given in the first row. 

"TP model or RHD snapshot {W) (LTE) (W) (NLTE) 



ID Kurucz (a = 0.5) 


38.4 


26.1 


T57G44N7-43108 


116.0 


27.7 


T57G44N7-69390 


115.4 


27.1 


T57G44N7-83752 


112.7 


27.3 


T57G44N7-87007 


110.9 


27.7 


T57G44N7-97644 


108.2 


26.5 



The mean equivalent widths derived from the horizontally averaged flux 
spectrum are listed in Table Q for the 5 snapshots and for the Kurucz ID reference 
model. In LTE, the 1D/2D difference is huge (granulation abundance correction 
~ —0.45 dex). But remarkably, the 2D NLTE line strengths show very little 
dispersion and do not indicate any significant offset with respect to the ID case. 
An obvious conclusion from these results is that the 2D LTE computations 



8 



R. Cayrel <fe M. Steffen 

OT57G44N7_69390x.f 1 5 



o 
o 



I I 

o o ^6 


^\ r 


O O o 

o o £*™saas 

. O o O 


§ o 
V <-7 ^ -~" ~x o ° / 




\ 7 \ <€> I 




o c VJ 

o — 




<# 




I I I I I 



o 



2 3 
x(Mm) 



WYLTE 
WYNLTE 
WYLTEmean = 
WYNLTEmean = 



95.6 



25.5 d£ o <g"" IO °e, 

o o o 

> ^ o CGDO CfSqp<iE> 

O <V 



o ' 

2x1 0" 



3x1 cr 



4x 1 0~ 



5x 1 0~ 



SURE(CGS) 



Figure 4. Same plots as Fig. 3, but for the snapshot from the metal- 
poor star simulation shown in Fig. 1 



are way off, strongly underestimating the Li abundance. In NLTE, the error 
introduced by representing the inhomogeneous stellar atmosphere by a flux- 
constant ID Kurucz model appears to be almost negligible. 

The mechanism behind the spatial variation of the line strength is clearly 
identified on Figs 1 and 4: hot granules produce at the same time a steeper 
temperature gradient and a lower temperature in the line formation region. In 
LTE, the latter leads to an overpopulation of the lower level of the transition 
due to a shift of the ionization equilibrium towards neutral particles (Saha equa- 
tion). Since both effects enhance the equivalent width of the line, the correlation 
between continuum intensity and LTE line strength is clearly positive (there is 
considerable dispersion around the mean W(I C ) relation, however, due to the 
presence of inclined thermal inhomogeneities) . This result is inverse of what is 
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expected for a set of hot and cool radiative equilibrium atmospheres lined up 
side by side: here the line would weaken in the hot atmosphere, because the 
population of the lower level is the dominant factor. It is clear that the reason- 
ing of Kurucz fails, essentially because the actual vertical stratification of hot 
and cold regions has little to do with the stratification of hot and cool radiative 
equilibrium atmospheres. 

A positive W-I c correlation alone is not sufficient to explain the LTE result 
that the mean equivalent width is much larger in 2D than in ID: as long as the 
fluctuations of the line strength with temperature remain in the linear domain, 
the mean equivalent width is not affected by the temperature fluctuations. But 
as the population of the ground level No of Li I varies exponentially with T, a 
nonlinearity sets in: symmetrical temperature fluctuations produce asymmet- 
rical population variations. Assume that Li is almost completely ionized and 
consider only the exponential factor of the Saha equation for simplicity. Then 

N (T) = N (T ) eM+x/kT) = N (T ) * jl - qs + qs 2 + X -q 2 s 2 + . . . J 

where x is the ionization potential and we have defined q = x/kTo and s = (T — 
To) /To, To being the mean temperature. Taking the horizontal average assuming 
symmetrical temperature variations ((s) = 0), we obtain (N) ~ No(To) * (1 + 
0.5q 2 {s 2 ) (y/ (s 2 ) < 0.04; q m 12). So the mean of W is biased towards larger 
equivalent widths in an inhomogeneous atmosphere. This explains part of the 
1D/2D difference found in our numerical computations. The main contribution, 
however, is attributed to the lower mean photospheric temperature in the 2D 
model as a result of adiabatic cooling due to overshooting: To (2D) < T(1D). 

In NLTE, the local temperature plays little role. Rather, the radiation field 
is the dominant factor. Hence, photoionization overionizes the lower level over 
a hot granule, and the equivalent width becomes smaller. As a result, NLTE 
equivalent widths are anti-correlated with the continuum surface brightness. 
The variation of W is much smaller than in LTE, because the angle-averaged 
radiation field depends only weakly on horizontal position at the height of line 
formation. 

Finally, we would like to mention the work by M. Asplund, who has also 
presented his NLTE line formation calculations for Li I on this symposium. 
His results are based on the 3D snapshots that he had used previously for the 
LTE investigation of the subgiants HD 140283 and HD 84937. He reached, on 
this completely independent set of hydro dynamical models, and with a different 
NLTE code, the same conclusions as we did: NLTE line formation in multidi- 
mensional models is quite different from the LTE case, in the way that the NLTE 
multi-dimensional Li abundance is much closer to the abundance derived from 
classical flux-constant ID models. A detailed comparison of our results shows 
that the remaining differences can be traced to the use of a Li I atom with 
5 levels in our case, and with 20 levels in the case of Asplund et al. (this vol- 
ume). Another difference is that the temperature inhomogeneities are somewhat 
enhanced in 2D models with respect of what occurs in 3D models, leading to 
correspondingly larger LTE granulation abundance corrections. But this is only 
a minor point. Certainly, the two groups agree that LTE Li abundance determi- 
nations relying on multidimensional hydrodynamical simulations of convection 
in metal-poor dwarfs are highly misleading. 
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6. Conclusions 

1. The statement of Kurucz (1995) that abundances of lithium derived from 
standard ID models of metal-poor stellar atmospheres is too small by a factor 
of 10 is not supported by actual multidimensional NLTE computations. Even 
the sign of the correction is doubtful, and the error is well below 0.1 dex, both 
according to our investigation and the one presented by M. Asplund on this 
symposium. 

2. LTE abundance determinations based on inhomogeneous atmospheres are 
strongly discouraged. They produce large "granulation abundance corrections" 
due to non-linear effects in the direction opposite to Kurucz' prediction, but the 
actual NLTE line formation mechanism couples the population of the atomic 
levels more closely to the mean radiation field than to the local temperature. 

3. The combination of multidimensional models with NLTE line formation for 
the Li I A 670.8 nm resonance line leads to the same lithium abundance as that 
derived from NLTE analysis with flux-constant ID models, abundance differ- 
ences being less than 0.1 dex. However, this result must not be hastily general- 
ized to other atoms with a different atomic structure. 

4. An obvious future improvement is to extend the NLTE analysis to the contin- 
uum, which has be assumed here to be in LTE. We plan to do that in the near 
future. If it turns out that the H~ ion is affected by NLTE, this will raise a new 
question: should such effects be included already in the radiation hydrodynamics 
code, which determines the amplitude of the thermodynamical fluctuations? 
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